home *** CD-ROM | disk | FTP | other *** search
/ The X-Philes (2nd Revision) / The X-Philes Number 1 (1995).iso / xphiles / hp48_1 / fftsrc_k < prev    next >
Internet Message Format  |  1995-03-31  |  21KB
  1. From: Ken Cooke <kcooke@milton.u.washington.edu>
  2. Subject:  v06i007:  fftsrc_kc - Machine Code FFT source v1.0, Part01/01
  3. Newsgroups: comp.sources.hp48
  4. Keywords: FFT
  5. Organization: University of Washington
  6. Followup-To: comp.sys.hp48
  7. Approved: spell@seq.uncwil.edu
  8.  
  9. Checksum: 1393729716 (verify with brik -cv)
  10. Submitted-by: Ken Cooke <kcooke@milton.u.washington.edu>
  11. Posting-number: Volume 6, Issue 7
  12. Archive-name: fftsrc_kc/part01
  13.  
  14.  
  15. BEGIN_DOC fftsrc.doc
  16. I'm glad to hear that someone liked the machine code FFT programs.
  17. In response to email requests, here is the SASM source code.  It can
  18. be assembled with the HP tools.  I didn't include FFT, because it is a
  19. subset of IFFT (I commented the lines to remove).
  20.  
  21. I apologize for the shortage of higher-level comments (most of this is still 
  22. on paper).  Also, I traded some readability to avoid subroutine calls (hardware
  23. stack levels are very precious) and keep variables in registers.  A few of
  24. the ROM calls are "unsupported entry points", which may need fixing in the
  25. future.  With rumors of new 48's in the works, ones with machine code Eqn
  26. Writers, re-written plotter user interfaces, and maybe even faster
  27. Saturn processors...well, lets just say that it may be wise to stick to
  28. entries listed in ENTRIES.A.
  29.  
  30. I hope this is of use to someone interested in ML math software.  While I
  31. like playing games on my 48 as much as the next guy, it would be nice to 
  32. see more math tools using the capabilities of sys-RPL and ML.  (I do, 
  33. occasionally, use my 48 as a calculating device :)  Unfortunately,
  34. there is no documentation on using the low-level floating point routines 
  35. (MULTF, RADDF, etc) so a little hacking with Voyager or MLDL is required.
  36.  
  37. I would love to hear from other math-hackers out there...
  38.  
  39. Ken Cooke
  40. N7VFE
  41. kcooke@u.washington.edu
  42. END_DOC
  43.  
  44.  
  45.  
  46. BEGIN_SRC angle.s
  47. * ANGLE.S  WORKS FOR ANY SIZE VECTOR OR ARRAY
  48. * Copyright 1992 Ken Cooke
  49. * Takes a VECTOR/ARRAY in STK1.  Returns a real VECTOR/ARRAY of
  50. * same size, that contains the angle of each complex element,
  51. * according to the current trig mode.
  52.  
  53.     TITLE ANGLE
  54. ASSEMBLE
  55.     NIBASC     /HPHP48-A/
  56. RPL
  57. ::
  58. CK1NoBlame
  59. CK&DISPATCH1
  60. FOUR
  61. ::
  62. * IF REAL ARRAY, EXIT UNCHANGED
  63.     DUP TYPECARRY? NOT?SEMI            ( ORIG )
  64. * MALLOC ANSWER ARRAY (SAME SIZE AS ORIG BUT REAL )
  65.     DUP MDIMS ITE TWO{}N ONE{}N %0 MAKEARRY    ( ORIG ANS )
  66. CODE
  67. ****************** UNLISTED ENTRY POINTS *****************
  68. ATTN?Lp    = #0CA88
  69. ARG15    = #2B6D7
  70. SetTrigMode    = #2AEF6
  71. ************************ MCODE ANGLE ***********************
  72. * ENTRY:
  73. *       STK2: ORIGINAL ARRY1 (COMPLEX)
  74. *    STK1: ANSWER ARRY (REAL)
  75.  
  76. * EXIT: STK1: ANSWER ARRAY (REAL)
  77.  
  78.     NIBHEX  823    ;CLEAR SB AND XM IN HST
  79.     D1=D1+    5    ;DROP STK1
  80.     D=D+1    A    ;RETURN MEM
  81.     GOSBVL    =SAVPTR    ;SAVE REGS WITH ONLY orig ON STACK
  82.  
  83.     A=DAT1    A    ;A->ORIG
  84.     D1=D1-    5    ;POINT BACK TO ANS
  85.     C=DAT1    A    ;C->ANS
  86.     R3=C.F    A    ;SAVE ->ANS IN R3
  87.     D0=C
  88.     D0=D0+    15
  89.     D0=D0+    10    ;D0->ANS DATA (OR DIM2 IF ARRAY)
  90.  
  91.     D1=A        ;D1->ORIG
  92.     D1=D1+    15    ;D1->#DIMS
  93.     A=DAT1    A    ;A=#DIMS
  94.     D1=D1+    5    ;D1->DIM1
  95.     C=DAT1    A    ;C=DIM1
  96.     D1=D1+    5    ;D1->ORIG DATA (OR DIM2 IF ARRAY)
  97.     A=A-1    A    ;CHECK FOR VECTOR
  98.     ?A=0    A
  99.     GOYES    ISVECT
  100. ***IF ARRAY, #ELEMENTS=DIM1*DIM2.  ALSO, SKIP OVER DIM2 FIELD
  101.     A=DAT1    A    ;A=DIM2
  102.     D1=D1+    5    ;D1->ORIG DATA
  103.     D0=D0+    5    ;D0->ANS DATA
  104.     SETHEX
  105.     GOSBVL    =MUL#    ;B=A*C
  106.     C=B    A    ;C=DIM1*DIM2
  107.  
  108. ISVECT  C=C-1    A    ;SUB 1 SO CAN USE BORROW FOR LOOPING
  109.     R4=C.F    A    ;SAVE #ELEMENTS-1 IN R4
  110.  
  111. ABSLOOP GOSBVL    =ATTN?Lp    ;ABORT IF ATTN (ON) WAS PRESSED
  112.     A=DAT1    W    ;A=RE
  113.     D1=D1+    16    ;D1->IM
  114.     CD0EX
  115.     RSTK=C        ;PUSH D0
  116.     C=DAT1    W    ;C=IM
  117.     D1=D1+    16    ;D1->NEXT RE
  118.     SETDEC
  119.     GOSBVL    =SPLTAC    ;CONVERT BOTH TO LONGS
  120.     GOSBVL    =SetTrigMode    ;SET ST ACCORDING TO TRIG MODE FLAGS
  121.     GOSBVL    =ARG15    ;AB=ANGLE (TRASHES A,B,C,D,D0,R0,R1)
  122.     GOSBVL    =PACK    ;A=ABSVAL
  123.     C=RSTK
  124.     D0=C        ;POP D0
  125.     DAT0=A    W    ;STORE ABSVAL IN ANSWER ARRAY
  126.     D0=D0+    16    ;POINT AT NEXT ANSWER
  127. * CHECK FOR DONE
  128.     C=R4.F    A    ;GET COUNTER
  129.     SETHEX
  130.     C=C-1    A
  131.     R4=C.F    A    ;DEC COUNTER
  132.     GONC    ABSLOOP    ;LOOP UNTIL BORROW
  133.  
  134. ***NOW ANSWER CONTAINS ABS VALUES OF ORIG
  135.     A=R3.F    A
  136.     GOSBVL    =GPOverWrALp    ;GETPTR, REPLACE TOS WITH ANS, LOOP
  137.  
  138. ENDCODE
  139. ;
  140. ;
  141. END_SRC
  142.  
  143.  
  144. BEGIN_SRC ifft.s
  145. ******************************************************************
  146. * MACHINE CODE IFFT by Ken Cooke
  147. * Copyright 1992
  148. ******************************************************************
  149.  
  150. * RADIX-2 DIT IFFT
  151. * Input: real or complex vector on STK1, with size=2^n
  152. * Ouput: complex vector containing IFFT
  153. * Note: uses the definition that FFT uses exp(-xxx) and IFFT uses exp(+xxx)
  154.  
  155.     TITLE IFFT
  156. ASSEMBLE
  157.     NIBASC     /HPHP48-A/
  158. RPL
  159. ::
  160. CK1NoBlame
  161. CK&DISPATCH1
  162. FOUR
  163. ::
  164.     DUP MDIMS             ( ARRY #N flag )
  165. * CHECK FOR MATRIX                  (NcaseDIMERR = 37DF6 )
  166.     NOT NcaseDIMERR               ( ARRY #N )
  167. * CHECK FOR SIZE=2^N
  168.     DUP                ( ARRY #N #N )
  169.  
  170. CODE
  171. ****************** CHECK THAT SIZE=POWER OF TWO **************
  172. *CHKSIZE2 REPLACES BINT (ON TOS) WITH TRUE IF BINT IS A
  173. * POWER OF 2, FALSE OTHERWISE.
  174. *NOTE: WORKS UP TO 2^19, AND RETURNS TRUE FOR BINT=0
  175. WrTLoop    EQU    #03B1A
  176. WrFLoop    EQU    #03B06
  177.     C=DAT1    A    ;C->BINT
  178.     CD1EX        ;SAVE D1 IN C.A
  179.     D1=D1+    5    ;D1->DATA
  180.     A=DAT1    A    ;A.A=BINT
  181.     D1=C        ;RESTORE D1
  182.  
  183.     SETHEX
  184.     C=A    A
  185.     A=-A-1    A    ;A=ONES COMP
  186.     C=-C    A    ;C=TWOS COMP
  187.     A=A!C    A    ;NO ZEROS IF BINT=2^N
  188.     A=A+1    A    ;WILL SET CARRY IF BINT=2^N
  189.     GOC    DoTrue
  190.     GOVLNG    WrFLoop    ;OVERWRITE BINT WITH FALSE AND CONT
  191. DoTrue    GOVLNG    WrTLoop    ;OVERWRITE BINT WITH TRUE AND CONT
  192. ENDCODE
  193.  
  194. *        ( ARRY #N FLAG )
  195. NcaseDIMERR    ( IF FALSE THEN "INVALID DIMENSION" ERROR )
  196.  
  197. * IF SIZE=1, RETURN WITH ARRY ON STACK
  198.     DUP#1= ITE DROP     ( ARRY #N )
  199.  
  200. ::
  201. * CREATE TEMP STORAGE (MUST BE 172 NIBS BIGGER THAN UNPACKED A FOR VARS)
  202.  
  203.     FORTYTWO OVER #* 172 #+ MAKEHEX$ SWAP    ( ARRY HEX$ #N )
  204.     ONE{}N C%0 MAKEARRY            ( ARRY HEX$ ANSARRY )
  205.  
  206. * DO THE FFT
  207.  
  208. CODE
  209. ****************** UNLISTED ENTRY POINTS *****************
  210. PUTAB1    = #2C066
  211. GETCD1    = #2C017
  212. aSIN15    = #2B6E0
  213. ATTN?Lp    = #0CA88
  214. longpi    = #2A45D
  215. varnibs    = 172
  216. ************************ MCODE FFT ***********************
  217. * ENTRY:
  218. *    STK3: (orig) real/complex data array
  219. *       STK2: (temp) hex$ (to hold vars and unpacked reals during FFT)
  220. *    STK1: (ans) complex answer array
  221.  
  222. * EXIT: STK1: complex FFT array (ans)
  223.  
  224.     NIBHEX  823    ;CLEAR SB AND XM IN HST
  225.     D1=D1+    10    ;DROP temp AND ans
  226.     D=D+1    A    ;RETURN MEM
  227.     D=D+1    A
  228.     GOSBVL    =SAVPTR    ;SAVE REGS WITH ONLY orig ON STACK
  229.  
  230. * SAVE NEEDED POINTERS IN SCRATCH REGS
  231.     D1=D1-    5    ;POINT BACK TO temp
  232.     A=DAT1    A    ;A->temp
  233.     A=A+CON    A,10    ;SKIP HEADER TO vars
  234.     R4=A.F    A    ;R4->vars
  235.     P=    0    ;*****REMEMBER:ALWAYS HAVE P=0 BEFORE LC(5)*****
  236.     LC(5)    varnibs    ;
  237.     SETHEX
  238.     A=A+C    A    ;SKIP OVER VARS
  239.     R1=A.F    A    ;R1->data
  240.  
  241.     D1=D1+    5    ;POINT TO orig
  242.     A=DAT1    A    ;A->orig
  243.     D0=A
  244.     D0=D0+    10    ;D0->OBJ PROLOG
  245.  
  246. * SET FLAG IF REAL, FOR USE WHEN UNPACKING
  247. * USE D.15 INSTEAD OF ST TO AVOID TROUBLE
  248.     D=0    S    ;USE D.15 AS FLAG (0=COMPLEX,1=REAL)
  249.     A=DAT0    A    ;A=OB PROLOG
  250.     LC(5)    #02977    ;C=COMPLEX PROLOG
  251.     ?A=C    A
  252.     GOYES    CMPLEX
  253.     D=D+1    S    ;IF REAL, D.15=1
  254.  
  255. CMPLEX    D0=D0+    10    ;D0->SIZE
  256.     A=DAT0    A    ;A.A=SIZE
  257.     R2=A.F    A    ;STORE SIZE IN R2
  258.     D0=D0+    5    ;D0->DATA
  259.     CD0EX
  260.     R0=C.F    A    ;R0->origdata
  261.  
  262. *FIND LOG2 OF SIZE
  263.     C=0    S    ;C.S WILL HOLD SIZE2=LOGbase2(SIZE)
  264. LOG2LP    ASRB.F    A
  265.     C=C+1    S
  266.     ?ABIT=0    0
  267.     GOYES    LOG2LP
  268.     R2=C.F    S    ;STORE SIZE2 IN R2.15
  269.  
  270. *************************** BITREV AND UNPACK ****************************
  271. *NOW, CONVERT REALS IN orig INTO EXT.REALS IN data
  272. * AT THE SAME TIME, REARRANGE INTO BIT-REVERSED ORDER
  273.  
  274.     A=R0.F    A
  275.     D0=A        ;D0->origdata
  276.     D=0    A    ;D HOLDS INDEX INTO orig (START AT 0)
  277.  
  278. BITREV    CPEX    15    ;PUT BIT-COUNTER IN P
  279.     P=P-1           ;SUB ONE SO CAN USE BORROW FOR BITLOOP
  280.     C=D    A    ;C GETS INDEX INTO orig
  281.     B=0    A    ;B WILL BE BIT-REV INDEX
  282.  
  283. *BIT REVERSAL
  284. BITLOOP B=B+B    A    ;SHIFT B LEFT ONE BIT
  285.     ?CBIT=0    0    ;TEST LSB
  286.     GOYES    NOCARRY    ;SKIP IF ZERO
  287.     B=B+1    A    ;ELSE FEED INTO B
  288. NOCARRY    CSRB.F    A       ;SHIFT C RIGHT ONE BIT
  289.     P=P-1        ;DEC BIT-COUNTER (BORROW WHEN DONE)
  290.     GONC    BITLOOP    ;LOOP IF MORE
  291.  
  292. *MULT B (BITREV INDEX) BY 21 AND ADD TO ->data
  293.     A=B    A    ;COULD USE =MUL# BUT THIS IS FASTER/SHORTER
  294.     B=B+B    A
  295.     B=B+B    A    ;B=4B
  296.     A=A+B    A    ;A=5B
  297.     B=B+B    A
  298.     B=B+B    A    ;B=16B
  299.     B=A+B    A    ;B=21B
  300.     B=B+B    A    ;B=42B
  301.  
  302.     A=R1.F    A    ;A->data
  303.     A=A+B    A    ;
  304.     D1=A           ;D1 POINTS INTO data AT BITREV POSITION
  305.     A=DAT0    W    ;GET RE(Ai)
  306.     D0=D0+    16    ;POINT AT NEXT REAL
  307.     GOSBVL    =SPLITA    ;NO NEED TO SETDEC (DONE IN SPLITA)
  308.     GOSBVL  =PUTAB1    ;PUT LONG AT DAT1, D1+=21
  309.  
  310. * IF INPUT ARRAY WAS REAL, USE IMAG=0
  311.     ?D=0    S    ;IF COMPLEX
  312.     GOYES    GETIM    ; GET IMAG PART
  313.     A=0    W    ;ELSE USE IMAG=0
  314.     B=0    W
  315.     GONC    NOGETIM    ;GO ALWAYS
  316.  
  317. GETIM    A=DAT0    W    ;GET IM(Ai)
  318.     D0=D0+    16
  319.     GOSBVL    =SPLITA
  320.  
  321. NOGETIM    GOSBVL  =PUTAB1    ;STORE IMAG PART
  322.     SETHEX
  323.     D=D+1    A    ;INC INDEX
  324.     C=R2        ;LOAD SIZE AND SIZE2 INTO C
  325.     ?D=C    A    ;COMPARE INDEX TO SIZE
  326.     GOYES    BITDONE
  327.     GOTO    BITREV    ;REPEAT IF NOT DONE
  328. BITDONE
  329.  
  330. ********** data NOW HOLDS LONG BITREV ARRAY **********
  331.  
  332.  
  333. ******************** BUTTERFLIES *******************
  334. * ENTRY: temp CONTAINS BIT-REV LONG COMPLEX DATA
  335. *        A IS NO LONGER NEEDED
  336. * Define template for var storage: (var names from Num.Recipes)
  337.  
  338. istep    = 0
  339. n    = (istep)+5
  340. mmax    = (n)+5
  341. data0    = (mmax)+5
  342. wr    = (data0)+5
  343. wi    = (wr)+21
  344. tempr    = (wi)+21
  345. wtemp    = (tempr)+21
  346. theta    = (wtemp)+21
  347. wpi    = (theta)+21
  348. wpr    = (wpi)+21
  349. m    = (wpr)+21
  350.  
  351. *FIRST, STORE SCRATCH REGS IN VARS
  352.  
  353.     SETHEX
  354.     A=R4.F    A    ;A->istep
  355.     D1=A
  356.     D1=D1+    5    ;D1->n
  357.     C=R2.F    A    ;C=size
  358.     C=C+C    A    ;C=2*size
  359.     DAT1=C    A    ;n=2*size
  360.  
  361.     D1=D1+    5    ;D1->mmax
  362.     P=    0
  363.     LC(5)    #2    ;C=2
  364.     DAT1=C    A    ;mmax=2
  365.  
  366.     D1=D1+    5    ;D1->data0    ***STORE POINTER TO data[0] HERE
  367.     C=R1.F    A    ;***THIS CORRECTS FOR "FORTRAN INDEXING"
  368.     C=C-CON    A,16    ;*** SINCE ALGORITHM ASSUMES INDEXES START AT 1
  369.     C=C-CON    A,5    ;*** I.E. data[1] = data0+21 = FIRST SPOT
  370.     DAT1=C    A
  371.  
  372. **REDUCED TO ONE SINE/MAINLOOP BY USING wtemp TO HOLD LAST SIN(theta)
  373. ** AND CALCULATING SIN(theta/2)
  374.     LC(5)    wtemp   ;INIT wtemp TO AVOID SIN(PI) SINCE NOT RETURN ZERO
  375.     C=C+A    A    ;A STILL =R4
  376.     D0=C        ;D0->wtemp
  377.     A=0    W
  378.     B=0    W    ;AB=0.0
  379.     GOSBVL    =PUTAB0    ;STORE wtemp=0.0, D0->theta
  380.  
  381. **INIT theta TO -pi FOR FFT, +pi FOR IFFT
  382.     D1=(5)  =longpi
  383.     GOSBVL    =GETAB1    ;get pi from ROM
  384. *    SETDEC                *********enable for FFT
  385. *    A=-A-1    S    ;A=-PI        *********enable for FFT
  386.     GOSBVL    =PUTAB0    ;STORE IN theta
  387.  
  388. ************************ BEGIN MAIN LOOP ***********************
  389. MAINLOOP
  390.     C=R4.F    A    ;C->vars
  391.     D1=C
  392.     D1=D1+    5    ;D1->n
  393.     C=DAT1    A    ;C=n
  394.     D1=D1+    5    ;POINT AT mmax
  395.     A=DAT1    A    ;A=mmax
  396. * MAINLOOP TEST
  397.     ?A<C    A    ;IF (mmax>=n)
  398.     GOYES    MAINCONT    ; THEN DONE
  399.     GOTO    MAINDONE
  400. ** A=mmax, D1->mmax
  401. MAINCONT
  402.     SETHEX
  403.     A=A+A    A    ;A=2*mmax
  404.     D1=D1-    10    ;istep
  405.     DAT1=A    A    ;STORE istep=2*mmax
  406.  
  407.     D1=D1+    10
  408.     D1=D1+    10    ;D1->wr
  409.     A=0    W
  410.     B=0    W
  411.     P=    14    ;CREATE AB=1.0
  412.     B=B+1    P
  413.     GOSBVL    =PUTAB1    ;STORE wr=1.0, D1->wi
  414.     B=0    P
  415.     GOSBVL    =PUTAB1    ;STORE wi=0.0, D1->tempr
  416.  
  417.     D1=D1+    16
  418.     D1=D1+    5    ;D1->wtemp
  419.     CD1EX
  420.     RSTK=C        ;PUSH ->wtemp
  421.     D0=C        ;D0->wtemp
  422.     D1=C        ;D1->wtemp
  423.     D1=D1+    16
  424.     D1=D1+    16
  425.     D1=D1+    10    ;D1->wpi
  426.  
  427.     SETDEC        :***NOTE: SETDEC BEFORE ALL FP MATH (WILL STAY IN DEC)
  428.     GOSBVL    =GETAB0    ;AB=wtemp, D0->theta
  429.     GOSBVL    =PUTAB1    ;STORE wpi=wtemp, D1->wpr
  430.     GOSBVL    =GETAB0    ;AB=theta, D0->wpi
  431.     GOSBVL    =DIV2    ;AB=theta/2
  432.     D0=D0-    16
  433.     D0=D0-    5    ;D0->theta
  434.     GOSBVL    =PUTAB0    ;STORE theta=theta/2
  435.     ST=1    9    ;SET UP FOR SIN (IN RADIANS)
  436.     ST=0    4
  437.     GOSBVL    =aSIN15    ;**NOTE:TRASHES A,B,C,D,P,D0,R0,R1!!!
  438.     C=RSTK        ;POP ->wtemp
  439.     D0=C        ;D0->wtemp
  440.     GOSBVL    =PUTAB0    ;STORE wtemp=SIN(theta)
  441.  
  442.     C=B    W
  443.     D=C    W    ;CD=AB
  444.     C=A    W
  445.     GOSBVL    =MULTF    ;SQUARE wtemp
  446.     C=B    W
  447.     D=C    W    ;CD=AB
  448.     C=A    W
  449.     GOSBVL    =RADDF    ;DOUBLE AB (D0 TRASHED)
  450.     A=-A-1    S    ;NEGATE AB
  451.     GOSBVL    =PUTAB1    ;STORE wpr=-2*SQR(wtemp), D1->m
  452.  
  453.     A=0    A
  454.     A=A+1    B    ;A=1
  455.     DAT1=A    A    ;STORE m=1
  456.  
  457. ************************ BEGIN OUTER LOOP ***********************
  458. OUTLOOP    SETHEX
  459.     A=R4.F    A
  460.     D0=A
  461.     D0=D0+    10    ;D0->mmax
  462.     P=    0
  463.     LC(5)    m
  464.     C=C+A    A
  465.     D1=C        ;D1->m
  466.     A=DAT0    A    ;A=mmax
  467.     C=DAT1    A    ;C=m
  468. *OUTLOOP TEST
  469.     ?C<A    A    ;IF (m>=mmax)
  470.     GOYES    NOTDONE    ; THEN DONE
  471.     GOTO    OUTDONE
  472. NOTDONE    RSTK=C        ;PUSH i=m (i KEPT ON TOS)
  473.     A=R4.F    A
  474.     D1=A        ;D1->n-5
  475.  
  476. ************************ BEGIN INNER LOOP ***********************
  477. INLOOP
  478.     GOSBVL    =ATTN?Lp    ;IF ATTN(ON) HAS BEEN PRESSED, GETPTR AND LOOP
  479.     D1=D1+    5    ;D1->n
  480.     A=DAT1    A    ;A=n
  481.     C=RSTK        ;C=i
  482. *INLOOP TEST
  483.     ?C<=A    A    ;IF (i>n)
  484.     GOYES    INCONT    ; THEN DONE
  485.     GOTO    INDONE
  486. INCONT    RSTK=C        ;PUSH i BACK
  487.     D1=D1+    5    ;D1->mmax
  488.     A=DAT1    A    ;A=mmax
  489.     SETHEX
  490.     A=A+C    A    ;A=j
  491.     B=A    A
  492.     B=B+B    A
  493.     B=B+B    A
  494.     A=A+B    A
  495.     B=B+B    A
  496.     B=B+B    A
  497.     A=A+B    A    ;A=21*j
  498.  
  499.     B=C    A
  500.     C=C+C    A
  501.     C=C+C    A
  502.     B=B+C    A
  503.     C=C+C    A
  504.     C=C+C    A
  505.     B=C+B    A    ;B=21*i
  506.  
  507.     D1=D1+    5    ;D1->data0
  508.     C=DAT1    A    ;C->data
  509.     A=A+C    A    ;A->data[j]
  510.     C=B+C    A    ;C->data[i]
  511.     RSTK=C        ;PUSH ->data[i]
  512.     C=A    A
  513.     RSTK=C        ;PUSH ->data[j]
  514.  
  515.     D1=D1+    5    ;D1->wr
  516.     D0=C        ;D0->data[j]
  517.     GOSBVL    =GETAB1    ;AB=wr, D1->wi
  518.     GOSBVL    =STAB2    ;SAVE wr IN R2/R3
  519.     GOSBVL    =GETCD0    ;CD=data[j], D0->data[j+1]
  520.     SETDEC
  521.     GOSBVL    =MULTF    ;AB=wr*data[j]
  522.     GOSBVL    =STAB0    ;SAVE IN R0/R1
  523.  
  524.     GOSBVL    =GETAB1    ;AB=wi, D1->tempr
  525.     GOSBVL    =GETCD0    ;CD=data[j+1]
  526.     GOSBVL    =MULTF    ;AB=wi*data[j+1]
  527.     GOSBVL    =RCCD0    ;CD=wr*data[j]
  528.     A=-A-1    S    ;NEGATE AB
  529.     GOSBVL    =RADDF    ;AB=RESULT (D0 TRASHED)
  530.  
  531.     GOSBVL    =PUTAB1    ;tempr=RESULT, D1->tempi
  532.     C=RSTK
  533.     RSTK=C        ;COPY ->data[j] FROM TOS
  534.     D0=C        ;D0->data[j]
  535.     GOSBVL    =GETCD0    ;CD=data[j], D0->data[j+1]
  536.     D1=D1-    16
  537.     D1=D1-    16
  538.     D1=D1-    10    ;D1->wi
  539.     GOSBVL    =GETAB1    ;AB=wi
  540.  
  541.     GOSBVL    =MULTF
  542.     GOSBVL    =STAB0    ;STORE wi*data[j] IN R0/R1
  543.     GOSBVL    =GETCD0    ;CD=data[j+1]
  544.     GOSBVL    =RCAB2    ;AB=wr
  545.     GOSBVL    =MULTF    ;AB=wr*data[j+1]
  546.     GOSBVL    =RCCD0    ;CD=wi*data[j]
  547.     GOSBVL    =RADDF    ;AB=tempi (D0 TRASHED)
  548.  
  549.     GOSBVL    =STAB0    ;SAVE tempi IN R0/R1
  550.     C=RSTK
  551. ***STORE D1=>data[j]
  552.     D1=C        ;D1->data[j]
  553.     C=RSTK
  554.     RSTK=C        ;COPY ->data[i] FROM TOS
  555.     C=C+CON    A,16
  556.     C=C+CON    A,5    ;C->data[i+1]
  557.     RSTK=C        ;PUSH ->data[i+1]
  558.     D0=C
  559.     GOSBVL    =GETCD0    ;CD=data[i+1]
  560.     GOSBVL    =STCD2    ;SAVE data[i+1] IN R2/R3
  561.  
  562.     A=-A-1    S    ;NEGATE tempi
  563.     GOSBVL    =RADDF    ;AB=data[i+1] - tempi
  564.     CD1EX
  565.     D1=C        ;C->data[j]
  566.     D0=C
  567.     D0=D0+    16
  568.     D0=D0+    5    ;D0->data[j+1]
  569.  
  570.     GOSBVL    =PUTAB0    ;STORE data[j+1]
  571.     GOSBVL    =RCAB0    ;AB=tempi
  572.     GOSBVL    =RCCD2    ;CD=data[i+1]
  573.     GOSBVL    =RADDF    ;AB=tempi+(data[i+1])
  574.     C=RSTK        ;POP ->data[i+1]
  575.     D0=C        ;D0->data[i+1]
  576.     GOSBVL    =PUTAB0    ;STORE data[i+1]
  577.  
  578.     A=R4.F    A
  579.     P=    0
  580.     LC(5)    tempr
  581.     SETHEX
  582.     A=A+C    A
  583.     D0=A        ;D0->tempr
  584.     GOSBVL    =GETAB0    ;AB=tempr
  585.     GOSBVL    =STAB0    ;STORE tempr IN R0/R1
  586.     C=RSTK
  587.     RSTK=C        ;COPY ->data[i] FROM TOS
  588.     D0=C        ;D0->data[i]
  589.     GOSBVL    =GETCD0    ;CD=data[i]
  590.     GOSBVL    =STCD2    ;SAVE data[i] IN R2/R3
  591.     SETDEC
  592.     A=-A-1    S    ;NEGATE tempr
  593.     GOSBVL    =RADDF    ;AB=data[i]-tempr
  594. **USE D1=>data[j]
  595.     GOSBVL    =PUTAB1    ;STORE data[j]
  596.  
  597.     GOSBVL    =RCAB0    ;AB=tempr
  598.     GOSBVL    =RCCD2    ;CD=data[i]
  599.     GOSBVL    =RADDF    ;AB=data[i]+tempr
  600.     C=RSTK        ;POP ->data[i]
  601.     D0=C        ;D0->data[i]
  602.     GOSBVL    =PUTAB0    ;STORE data[i]
  603.  
  604.     A=R4.F    A
  605.     D1=A        ;D1->istep
  606.     A=DAT1    A    ;A=istep
  607.     C=RSTK        ;POP i
  608.     SETHEX
  609.     C=C+A    A    ;i=i+istep
  610.     RSTK=C        ;PUSH NEW i
  611.  
  612.     GOTO    INLOOP    ;NEXT ITERATION
  613. ************************ END INNER LOOP ***********************
  614. INDONE
  615.     A=R4.F    A
  616.     P=    0
  617.     LC(5)    wr
  618.     SETHEX
  619.     C=C+A    A
  620.     RSTK=C        ;PUSH *wr
  621.     D0=C        ;D0->wr
  622.     LC(5)    wpr
  623.     C=C+A    A
  624.     D1=C        ;D1->wpr
  625.  
  626.     SETDEC
  627.     GOSBVL    =GETAB0    ;AB=wr
  628.     GOSBVL    =STAB0    ;SAVE OLDwr IN R0/R1
  629.     GOSBVL    =GETCD1    ;CD=wpr, D1->wpr+21
  630.     GOSBVL    =MULTF    ;AB=wr*wpr
  631.     GOSBVL    =STAB2    ;SAVE IN R2/R3
  632.     GOSBVL    =GETCD0    ;CD=wi
  633.     D1=D1-    16
  634.     D1=D1-    16
  635.     D1=D1-    10    ;D1->wpi
  636.     GOSBVL    =GETAB1    ;AB=wpi, D1->wpr
  637.     GOSBVL    =MULTF    ;AB=wi*wpi
  638.     A=-A-1    S    ;NEGATE AB
  639.     GOSBVL    =RCCD2    ;CD=wr*wpr
  640.     GOSBVL    =RADDF    ;AB=wr*wpr - wi*wpi (D0 TRASHED)
  641.     GOSBVL    =RCCD0    ;CD=OLDwr
  642.     GOSBVL    =RADDF    ;AB=RESULT (D0 TRASHED)
  643.  
  644.     C=RSTK        ;POP *wr
  645.     D0=C        ;D0->wr
  646.     GOSBVL    =PUTAB0    ;STORE NEW wr, D0->wi
  647.     D1=D1-    16
  648.     D1=D1-    5    ;D1->wpi
  649.     GOSBVL    =GETAB1    ;AB=wpi, D1->wpr
  650.     GOSBVL    =RCCD0    ;CD=OLDwr
  651.     GOSBVL    =MULTF    ;AB=OLDwr*wpi
  652.     GOSBVL    =STAB0    ;SAVE IN R0/R1
  653.  
  654.     GOSBVL    =GETCD0    ;CD=wi
  655.     GOSBVL    =STCD2    ;SAVE wi IN R2/R3
  656.     GOSBVL    =GETAB1    ;AB=wpr, D1->m
  657.     GOSBVL    =MULTF    ;AB=wpr*wi
  658.     GOSBVL    =RCCD0    ;CD=OLDwr*wpi
  659.     GOSBVL    =RADDF    ;AB=wpr*wi + OLDwr*wpi (D0 TRASHED)
  660.     GOSBVL    =RCCD2    ;CD=wi
  661.     GOSBVL    =RADDF    ;AB=RESULT
  662.  
  663.     C=R4.F    A
  664.     D0=C
  665.     D0=D0+    16
  666.     D0=D0+    16
  667.     D0=D0+    9    ;D0->wi
  668.     GOSBVL    =PUTAB0    ;STORE NEW wi
  669.     C=DAT1    A    ;C=m
  670.     SETHEX
  671.     C=C+1    A
  672.     C=C+1    A    ;m=m+2
  673.     DAT1=C    A    ;STORE NEW m
  674.  
  675.     GOTO    OUTLOOP    ;NEXT ITERATION
  676. ************************ END OUTER LOOP ***********************
  677.  
  678. OUTDONE    SETHEX        ;D0 STILL POINTS TO mmax, A=mmax
  679.     A=A+A    A    ;A=2*mmax
  680.     DAT0=A    A    ;STORE NEW mmax
  681.  
  682.     GOTO    MAINLOOP
  683. ************************ END MAIN LOOP ***********************
  684.  
  685. MAINDONE
  686. * temp NOW HOLDS vars (FIRST 172 NIBS) AND LONG COMPLEX ARRAY OF ANSWERS
  687.  
  688. ***************** PACK AND COPY INTO ans *****************
  689.  
  690.     A=R4.F    A
  691.     D1=A        ;D1->vars
  692.     D1=D1+    5    ;D1->n
  693.     A=DAT1    A    ;A=n
  694.     R2=A.F    A    ;SAVE n IN R2 (COUNTER IN PACKLOOP)
  695.  
  696. ***REMOVE FOR FFT
  697.     ASRB.F    A    ;A=size
  698. ************* CONVERT A.A INTO EXT.REAL ************
  699. * USING HORNER'S RULE, 10100 = 2^4 + 2^2 = ((((1)2+0)2+1)2+0)2+0
  700. *  SINCE THE *2 AND ADDITION IS DONE IN DEC, RESULT IN DEC
  701. * ENTER: BINT IN A.A (ANY)
  702. * EXIT: EXT.REAL IN A/B (DEC)
  703. * TRASHES: A,B,C,P
  704.  
  705.     P=    0
  706.     LCHEX  19    ;REPEAT FOR EACH BIT IN A.A
  707.     B=0    W    ;B WILL HOLD MANTISSA
  708.  
  709. BITLP    B=B+B  W    ;B*=2
  710.     SETHEX
  711.     A=A+A  A    ;SHIFT MSB OF BINT INTO CARRY
  712.     SETDEC
  713.     GONC   SKIP
  714.     B=B+1  W
  715. SKIP    C=C-1  B    ;DEC BIT COUNTER
  716.     GONC   BITLP
  717.  
  718.     P=     5    ;A.A IS 0 FROM SHIFTING IN ZEROS
  719.     A=A-1  X
  720. NORMLP    A=A+1  X    ;A ACCUM'S EXPONENT BY COUNTING DIGITS
  721.     BSRC        ;ROTATE UNTIL MANT IS LEFT-JUSTIFIED
  722.     ?B#0   WP
  723.     GOYES  NORMLP
  724.  
  725.     BSR    W    ;SHIFT MANT INTO PLACE (B.S=0)
  726.     A=0    S    ;SIGN IS POSITIVE
  727. * AB NOW HOLDS EXT.REAL
  728. **************** END OF CONVERSION ****************
  729.     GOSBVL    =STAB0    ;SAVE SIZE IN R0/R1
  730. *******REMOVE FOR FFT
  731.  
  732.     SETHEX
  733.     A=R4.F    A    ;**THIS IS THE LAST TIME WE USE R4(temp)
  734.         P=    0
  735.     LC(5)    varnibs    ;SKIP OVER VARS TO DATA
  736.     C=A+C    A    ;C->data
  737.     D0=C        ;D0->data (SOURCE OF COPY)
  738.  
  739.     D1=(5)    =DSKTOP    ;RECALL SAVED D1
  740.     C=DAT1    A
  741.     D1=C        ;D1->STK3 (orig)
  742.     D1=D1-    10    ;D1->STK1 (ans)
  743.     C=DAT1    A
  744.     R3=C.F    A    ;SAVE ->ans IN R3
  745.     D1=C        ;D1->ans
  746.     D1=D1+    15      ;SKIP HEADER
  747.     D1=D1+    10    ;D1->ansdata (DEST OF COPY)
  748.  
  749. * NOW PACK AND COPY
  750. * NOTE: OVER/UNDERFLOW ERRORS ARE TRAPPED IN =PACK. orig IS SAVED AS TOS
  751. PAKLOOP
  752.     GOSBVL    =GETAB0    ;GET NEXT LONG
  753.     GOSBVL    =RCCD0    ;RECALL SIZE        *REMOVE FOR FFT
  754.     SETDEC
  755.     GOSBVL    =DIVF    ;DIVIDE BY SIZE        *REMOVE FOR FFT
  756.     GOSBVL    =PACK    ;EXT.REAL->REAL
  757.  
  758.     DAT1=A    W    ;PUT REAL IN PLACE
  759.     D1=D1+    16    ;POINT AT NEXT REAL
  760.  
  761.     C=R2.F    A    ;GET COUNTER
  762.     SETHEX
  763.     C=C-1    A    ;DEC COUNTER
  764.     R2=C.F    A
  765.     ?C#0    A
  766.     GOYES    PAKLOOP
  767.  
  768.  
  769. * NOW RESULT IS PACKED INTO ans ARRAY
  770. * AT LAST, REPLACE STK1 (ORIGINAL ARRAY) WITH RESULT ARRAY
  771.     A=R3.F    A        ;A->ans
  772.     GOVLNG    =GPOverWrALp    ;RESTORE REGS, PUT @A ON TOS, CONT RPL
  773. ENDCODE
  774. ;
  775. ;
  776. ;
  777. END_SRC
  778.  
  779.  
  780.  
  781. BEGIN_SRC absv.s
  782. * ABSV.S  WORKS FOR ANY SIZE VECTOR OR ARRAY
  783. * Copyright 1992 Ken Cooke
  784. * Takes a VECTOR/ARRAY in STK1.  Returns a real VECTOR/ARRAY of
  785. * same size, that contains the element-by-element absolute values.
  786.  
  787.     TITLE ABSV
  788. ASSEMBLE
  789.     NIBASC     /HPHP48-A/
  790. RPL
  791. ::
  792. CK1NoBlame
  793. CK&DISPATCH1
  794. FOUR
  795. ::
  796. * IF REAL ARRAY, EXIT UNCHANGED
  797.     DUP TYPECARRY? NOT?SEMI            ( ORIG )
  798. * MALLOC ANSWER ARRAY (SAME SIZE AS ORIG BUT REAL )
  799.     DUP MDIMS ITE TWO{}N ONE{}N %0 MAKEARRY    ( ORIG ANS )
  800. CODE
  801. ****************** UNLISTED ENTRY POINTS *****************
  802. ATTN?Lp    = #0CA88
  803. SQR15    = #2B9F3
  804.  
  805. ************************ MCODE ABSV ***********************
  806. * ENTRY:
  807. *       STK2: ORIGINAL ARRY1 (COMPLEX)
  808. *    STK1: ANSWER ARRY (REAL)
  809.  
  810. * EXIT: STK1: ANSWER ARRAY (REAL)
  811.  
  812.     NIBHEX  823    ;CLEAR SB AND XM IN HST
  813.     D1=D1+    5    ;DROP STK1
  814.     D=D+1    A    ;RETURN MEM
  815.     GOSBVL    =SAVPTR    ;SAVE REGS WITH ONLY orig ON STACK
  816.  
  817. * NOW IF USER HITS "ON" KEY, WILL LEAVE ORIGINAL STACK
  818.     A=DAT1    A    ;A->ORIG
  819.     D1=D1-    5    ;POINT BACK TO ANS
  820.     C=DAT1    A    ;C->ANS
  821.     R0=C.F    A    ;SAVE ->ANS IN R0
  822.     D0=C
  823.     D0=D0+    15
  824.     D0=D0+    10    ;D0->ANS DATA (OR DIM2 IF ARRAY)
  825.  
  826.     D1=A        ;D1->ORIG
  827.     D1=D1+    15    ;D1->#DIMS
  828.     A=DAT1    A    ;A=#DIMS
  829.     D1=D1+    5    ;D1->DIM1
  830.     C=DAT1    A    ;C=DIM1
  831.     D1=D1+    5    ;D1->ORIG DATA (OR DIM2 IF ARRAY)
  832.     A=A-1    A    ;CHECK FOR VECTOR
  833.     ?A=0    A
  834.     GOYES    ISVECT
  835. ***IF ARRAY, #ELEMENTS=DIM1*DIM2.  ALSO, SKIP OVER DIM2 FIELD
  836.     A=DAT1    A    ;A=DIM2
  837.     D1=D1+    5    ;D1->ORIG DATA
  838.     D0=D0+    5    ;D0->ANS DATA
  839.     SETHEX
  840.     GOSBVL    =MUL#    ;B=A*C
  841.     C=B    A    ;C=DIM1*DIM2
  842.  
  843. ISVECT  C=C-1    A    ;SUB 1 SO CAN USE BORROW FOR LOOPING
  844.     R4=C.F    A    ;SAVE #ELEMENTS-1 IN R4
  845.  
  846. ABSLOOP GOSBVL    =ATTN?Lp    ;ABORT IF ATTN (ON) WAS PRESSED
  847.     SETDEC
  848.     GOSUB    SQDAT1    ;SQUARE RE PART
  849.     GOSBVL    =STAB2    ;SAVE IN R2/R3
  850.     GOSUB    SQDAT1    ;SQUARE IM PART
  851.  
  852.     CD0EX
  853.     RSTK=C        ;PUSH D0
  854.     GOSBVL    =RCCD2    ;CD=SQR(RE)
  855.     GOSBVL    =RADDF    ;AB=SQR(RE)+SQR(IM) (D0 TRASHED)
  856.     GOSBVL    =SQR15    ;AB=LONG ABSVAL
  857.     GOSBVL    =PACK    ;A=ABSVAL
  858.     C=RSTK
  859.     D0=C        ;POP D0
  860.     DAT0=A    W    ;STORE ABSVAL IN ANSWER ARRAY
  861.     D0=D0+    16    ;POINT AT NEXT ANSWER
  862. * CHECK FOR DONE
  863.     C=R4.F    A    ;GET COUNTER
  864.     SETHEX
  865.     C=C-1    A
  866.     R4=C.F    A    ;DEC COUNTER
  867.     GONC    ABSLOOP    ;LOOP UNTIL BORROW
  868.  
  869. ***NOW ANSWER CONTAINS ABS VALUES OF ORIG
  870.     GOVLNG    =GPOverWrR0Lp    ;GETPTR, REPLACE TOS WITH ANS, LOOP
  871.  
  872. ************************ SQDAT1 ***********************************
  873. * GET REAL FROM DAT1, UNPACK AND SQUARE IT, INC D1
  874.  
  875. SQDAT1    A=DAT1    W    ;A=REAL
  876.     D1=D1+    16    ;D1->NEXT
  877.     GOSBVL    =SPLITA    ;AB=REAL
  878.     C=B    W
  879.     D=C    W
  880.     C=A    W    ;COPY INTO CD
  881.     GOSBVL    =MULTF    ;AB=SQR(REAL)
  882.     RTN
  883.  
  884. ENDCODE
  885. ;
  886. ;
  887. END_SRC
  888.